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ABSTRACT 


The  present  paper  provides  « methodology  for  the  computation  of  the 
hlt/klll  probability  of  a high  caliber  weapon  shooting  In  a forestA 
We  assume  that  trees  are  distributed  at  random  according  to  a Poisson 
Law.  Bullets  which  hit  trees  with  sufficiently  high  velocity  can 
penetrate  the  trunks  and  continue  towards  the  mrget.  A model  Is 
provided  for  the  determination  of  the  distribution  of  the  exit  velocity 
of  a bullet.  Recursive  method  Is  given  for  the  computation  of  succes- 
sive exit  distributions  as  functions  of  the  Initial  (muzzle)  velocity, 
the  distances  between  the  trees  and  their  characteristics.  On  the 
basis  of  this  recursive  method  the  kill  probabilities  are  computed. 
Numerical  examples  are  provided  as  well  as  FORTRAN  programs. 
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Shooting  In  Forests,  Exit  Velocity,  Distributions,  Poisson  distribution, 
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The  present  technical  report  provides  a method  of  evaluating  the  degrada- 
tion effects  on  large  caliber  weapons  located  In  a forest  and  shooting 
at  a target  In  the  forest.  The  degradation  effect  is  actually  measured 
by  the  decrease  In  the  hit  and  kill  probabilities  of  these  weapons, 
compared  to  those  when  there  are  no  obstacles  In  the  trajectories  of 
the  bullets.  The  obstacles  considered  here  are  randomly  located  trees 
of  varying  trunk  size.  The  Initial  (muzzle)  velocity  of  the  bullets  is 
sufficiently  high  to  allow  penetration  through  trees.  However,  the 
penetrating  bullet  loses  energy  and  Its  exit  velocity  may  be  considerably 
smaller  than  the  velocity  at  which  it  hits  the  tree.  The  hitting  velocity 
depends  on  the  Initial  velocity  and  on  the  distance  of  the  tree  from  the 
origin.  The  exit  velocity  Is,  however,  a random  variable  which  depends 
on  the  length  of  the  bullet's  path  through  the  trunk  and  the  resistance 
of  the  wood  (type  of  tree).  The  penetrating  bullet  may  also  be  deflected 
from  Its  original  aimed  path.  In  Section  2 we  specify  the  physical  model 
under  consideration  and  the  probabilistic  assumptions.  The  distribution 
of  the  exist  velocity  Is  derived  from  this  model  analytically  in  Section  3. 
In  Section  4 we  develop  a recursive  method  for  the  numerical  determination 
of  the  consecutive  distributions  of  the  exist  velocities  from  n trees 
(n>l),  which  are  located  at  specified  distances  dj,  dg,...^  from  each 
other,  on  the  path  of  the  bullet.  This  recursive  method  Is  particularly 
convenient  for  numerical  analysis.  In  Section  5 we  derive  an  analytic 
expression  for  the  hit/kill  probability  of  a round.  This  analytic 
expression  requires,  however,  numerical  methods  of  evaluation.  We 
provide  numerical  methods  which  combine  exact  computations  with  some 
Monte  Carlo  estimation.  This  Monte  Carlo  estimation  appears  only  In  one 
stage  of  the  numerical  evaluation,  replacing  a complicated  numerical 
Integration.  The  method  Is  not,  however,  a pure  simulation  procedure. 

It  has  the  property  that  with  a very  small  number  of  (Independent)  runs 
we  attain  estimates  with  very  high  precision.  Two  alternative  procedures 
are  compared  with  respect  to  their  precision  and  required  computing 
time.  FORTRAN  programs  are  provided  In  the  appendices. 
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2.  THE  PHYSICAL  AND  THE  PROBABILISTIC  MODEL 

Consider  a weapon  located  at  the  origin,  0,  and  shooting  at  a target 
which  Is  at  range  R [m].  The  Initial  (muzzle)  velocity  of  the  bullet  Is 
vQ  [m/sec].  If  there  are  no  obstacles  along  the  trajectory  of  the 
bullet.  It  will  hit  the  target  with  probability  PH(vQ,R).  Given 
that  the  bullet  hits  the  target,  the  kill  probability  depends  on  the 
velocity  of  hitting  the  target,  which  Is  vqx(v0,R)  and  on  other  possible 
factors.  Let  PK(vQ,R)  designate  the  combined  kill  probability.  This 
function  is  specified  In  each  particular  case  according  to  the  specific 
weapon  and  fighting  conditions.  Similarly  the  function  x(vQ,R)  depends 
on  the  type  of  weapon,  etc.  We  will  consider  here,  for  the  sake  of 
simplicity,  a linear  decreasing  function  of  R,  vQ  x(vQ,R)  » vq-bR. 

This  Is  a good  approximation  when  the  Initial  velocity,  vQ,  Is  high  and 
R Is  not  too  large  a fraction  of  the  weapons  maximum  range.  The  method 
developed  In  the  present  paper  can  be  easily  generalized  to  other  types 
of  ballistic  functions. 

The  problem  of  shooting  In  the  forest  Is  that  of  randomly  placed  obstacles 
(trees)  along  the  path  (trajectory)  of  the  bullet.  We  assume  that  the 
trees  are  randomly  located  according  to  a Poisson  Law,  with  a given 
density,  p[No.  of  trees/m^].  Thus,  If  we  consider  a strip  around  the 
straight  line  connecting  the  origin  with  the  target,  of  length  R[m] 
and  width  l[m],  the  number  of  trees  to  be  found  on  this  strip  is  a random 
variable,  N,  having  a Poisson  distribution  with  mean  pR. 

In  case  of  N * n,  n^,  1,  let  Dj,  D2 Dn  be  the  distances  (In  [m])  from 

the  origin  to  the  location  of  the  center  of  the  first  tree;  from  the  first 
tree  to  the  center  of  the  second,  etc. 
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It  Is  well  known  that  the  location  points  (Figure  1)  { ■ D^,  • Dj  ♦ 0o> 

....  ■ Dj  ♦ ...  ♦ Dn,  are  random  variables  having  a joint  distribution 

like  the  order  statistics  In  a sample  of  n Independent  and  Identically 
distributed  random  variables  from  a uniform  distribution  on  [O.R] 

(H.  A.  David  [1]  pp.  80). 


The  trees  are  generally  of  varying  size.  We  are  actually  concerned  with 
the  size  of  the  trunk  at  a certain  height,  h,  above  the  ground.  For  the 
purpose  of  modeling  we  assume  that  a cut  along  a horizontal  plan  yields 
a circle  of  radius  T (Figure  2).  This  radius  Is  generally  a random 
variable  with  a specific  distribution,  Fy(t).  The  methods  developed 
In  the  present  paper  are  for  a fixed  radius  T.  We  discuss  at  the  end 
how  the  numerical  procedures  can  be  extended  to  cover  the  case  of  varying 
radius.  Notice  that  the  length  of  the  bullets  path  In  the  trunk  Is 


2 (T2  - u2)1/2 


(2.1) 


where  u (the  distance  of  the  path  from  the  center)  Is  a random  variable, 
having  a uniform  distribution  on  the  Interval  (0,T).  Accordingly,  l Is 
a random  variable  having,  for  a given  T,  a distribution  function  (c.d.f.) 


(x;T) 


' [‘  * &] 


0 < x < 2T. 


(2.2) 


3.  THE  DISTRIBUTION  OF  THE  EXIT  VELOCITY 


Given  an  Initial  velocity  vQ  and  a tree  of  radius  T located  at  a distance 


D from  the  origin,  the  question  Is  what  Is  the  distribution  of  the  exit 
velocity,  Vj,  of  a bullet  going  through  that  tree.  We  say  that  the 
exit  velocity  Is  zero  If  the  bullet  Is  absorbed  In  the  tree. 
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The  basic  physical  equation  Is  that  of  energy  conservation,  namely 


m 
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(\(v0.d)v0)2 


♦ u. 


(3.1) 


where  d ■ D - T,  m is  the  mass  of  the  bullet,  vQ  MvQ,d)  Is  the  entrance 
velocity,  l the  length  of  the  bullet's  path  within  the  tree  and  y a proper^ 
constant,  which  depends  on  the  tree's  resistance  (In  units  of  [g][m]/[sec]  ). 
The  physical  model  (3.1)  Is  accepted  to  be  a good  first  approximation  to 
the  more  complicated  phenomenon  of  a projectile  penetrating  a solid  mass 
(see  Lambert  [3]).  From  (3.1)  we  can  write 


vl2  ’ <vo  *<Vd»2  * aL’  (3*2) 

p 

where  a Is  a proper  constant.  Let  L*(vQ,  d)  • (V0X(V0»  /<x'  If  1 -L* 

the  bullet  will  be  absorbed  In  the  tree.  The  probability  of  this  event 
Is,  according  to  (2.2) 


q(vQ,d;  T) 


0 

. If  L*(vQ,d)  > 2T 

[, 

1/2 

(3.3) 

1 

, otherwise 

Similarly,  the  c.d.f.  of  the  exit  velocty,  Vj,  Is 


i1  V T) 


1 

1 - 


((*„  \(v„.d))2  - »!2)J 


0 " O’ 
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V1  > v0  x(vo*  d) 
o ivl  < vo  *(Vd) 


(3.4) 
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where  C*3+  ■ max  (a,o).  Notice  that  If  T Is  small,  vQ  x (v  ,d)  - v1 
may  be  greater  than  2aT. 

In  these  cases  the  c.d.f.  value  Is  zero.  Thus,  let 


'l#(Vd;T)  "[ 


(V0X(  V0*^)  ) ■ 2aT 


11/2 
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If  q(vQ,d;T)  ■ 0 then  (vQ,d;T)  .>  0°  the  other  hand.  If  q(vQ,d;T)  > 0 

then  v1*(v(),d;T)  ■ 0.  These  relationships  are  Illustrated  In  Figure  3. 

4.  RECURSIVE  DETERMINATION  OF  THE  SUCCESSIVE  EXIT  VELOCITY  DISTRIBUTIONS 

In  the  following  we  will  adopt  the  simple  model  v0MvQ,d)  ■ vQ-fld. 

This  assumption  Is  not  restrictive.  The  following  formulae  are  developed 
for  this  particular  function,  since  the  data  showed  linearity  In  the  region 
of  Interest.  The  formulae  can  be  easily  modified  for  other  types  of 
ballistic  functions.  Define,  Hj(x;vQ,d,T)  ■ H(x;vQ,d,T),  0^  x^  v-fid. 

4.1  The  Case  of  n»2 

Let  dj  and  d2  be  the  given  distances.  The  distribution  of  the  exit 
velocity  from  the  second  tree,  V2,  can  be  obtained  from  H^xjVQ.d.T) , 
since  the  exit  velocity  from  the  first  tree,  If  given,  can  be  applied 
to  compute  the  entrance  velocity  Into  the  second  tree.  Accordingly,  the 
c.d.f.  of  Vg,  given  vQ,  d^  d2  and  T Is 


v*? 


h2(v;  Vdi»d2*T) 


J Hi(v»x»d2» 


T)dH1(x;v0d1,T) 


(4.1) 
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Thus,  we  obtain  after  some  manipulations 


H2*v; Vdl,d2,T)  " 


[(v0-Bdred2)2  - v2]2 


• 

2«2T2  v 


V6?!  -Bd2 


[(v^Bdj)2  - (y+8d2)2]2 

1 U? 


1/2 


1/2 


(4.2) 


l (y2  - v2)2"!  "2 

I??- 


This  Integral  should  be  understood  as  a regular  Integral  over  the  range 
of  values  over  which  the  functions  within  the  squared  brackets,  [ ],  are 
both  positive.  Furthermore, 


(0  , if  G(y)  < 0 

<?(y)dy  , If  G(y)  > 0, 

where  G*(y)  is  the  derivative  of  G(y).  Our  approach  is  to  evaluate  (4.2) 
numerically.  We  therefore  leave  it  in  Its  present  form,  without  further 
analytical  elaboration.  For  the  purpose  of  approximating  H2(v;vQ,d1 ,d2,T) 
numerically  we  partition  the  Interval  (v,  v^Bd^-Bdg)  Into  M subintervals 
of  equal  size  A ■ (vo-Bd^-Bd2-v)/M. 
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As  M grows  (to  infinity)  the  right  hand  side  of  (4.4)  approaches  that  of 
(4.2). 

In  Table  1 we  present  the  results  of  some  computations  of  the  c.d.f. 

H2(v;Vo,di,d2,T)  according  to  approximation  (4.4).  These  computations 
were  performed  according  to  Program  BULLET  of  Appendix  1,  with  the  proper 
parameters  and  M * 400  and  M = 500.  We  have  tried  the  approximation  also 
with  M = 50,  but  for  small  T values  (0.1  and  0.2)  it  has  not  yielded 
accurate  results  for  small  v values. 

4.2  The  General  Case 

After  computing  the  values  of  H2(v;v0,d1,d2,T) , at  specified  values  of  v, 
one  can  compute  H3(v;  v0,d1,d2,d3,T)  at  those  values  of  v,  etc.  The 
computation  is  based  on  the  recursive  formula 


H 


n 


H1(v;x,d[1,T)  dHn_1(x;v0,d(n'1),T) 


(4.5) 


n-1 

v0-e  £ dj 

3 Hn-l(v+Bdn’V£(n~1)j)  + /J=1  ^(vix^.T)  dH^  (x;v(),d(n“1  },T) 

V+6<L 


n 

For  every  o _<  v _<  v -e  £ d . , 
0 j-1  J 


where  d(n“1}  - (dx d^  d(n)  = (i(n_1),dn) 
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TABLE  1.  The  distributions  H2(v;  vQ,  d2,  T)  for  vQ  • 1300  [m/sec], 
0 ■ 1.8,  a ■ 2,812,500  [m/sec2],  dx  - 200,  d2  ■ 250. 


V/T 

.1 

.2 

.3 

.4 

.5 

0. 

0. 000000 

0.743070 

0.966422 

0. 990392 

0.996216 

100. 

0. 000000 

0.743149 

0.967444 

0.990663 

0.996321 

200. 

0. 000000 

0.730490 

0.970347 

0.991447 

0.996623 

300. 

0.000000 

0.731469 

0.974683 

0.992632 

0.997082 

400. 

0. 000000 

0.813191 

0.979840 

0.994068 

0.997642 

300. 

0. 000000 

0.364330 

0.983136 

0.993378 

0.998234 

600. 

0. 000000 

0.363422 

0.990042 

0.996997 

0.998793 

700. 

0.237341 

0.963970 

0.994063 

0.998188 

0.999269 

300. 

0.373220 

0.983473 

0.996987 

0.999070 

0.999623 

900. 

0.792740 

0.993628 

0.998799 

0.999626 

0.999848 

1000. 

0.969721 

0.998369 

0.999683 

0.999901 

0.999960 

1100. 

0.997259 

0.999833 

0.999968 

0.999990 

0.999996 

1200. 

0.999998 

1. 000000 

1.000000 

1. 000000 

1.000000 

PROCESSOR  USAGE*  102.4  UNITS 


V/T 

.1 

.2 

.3 

.4 

.5 

Q 

m ■ cnn 

0.000000 

0.713346 

0.966422 

0.990392 

0.996216 

M 500  100. 

0. 000000 

0.716970 

0.967443 

0.990663 

0.996321 

200. 

0. 000000 

0.731379 

0.970347 

0.991447 

0.996623 

300. 

0. 000000 

0.763288 

0.974683 

0.992632 

0.997032 

400. 

0. 000000 

0.796374 

0.979840 

0.994068 

0.997642 

300. 

0. 000000 

0.841691 

0.983136 

0.993373 

0.998234 

600. 

0. 000000 

0.881748 

0.990042 

0.996997 

0.998793 

700. 

0.239113 

0.963970 

0.994063 

0.998183 

0.999269 

800. 

0.336019 

0.983473 

0.996987 

0.999070 

0.999623 

900. 

0.794083 

0.993628 

0.998799 

0.999626 

0.999848 

1000. 

0.969721 

0.998369 

0.999683 

0.999901 

0.999960 

1100. 

0.997239 

0.999833 

0.999963 

0.999990 

0.999996 

1200. 

0.999998 

1.000000 

1. 000000 

1. 000000 

1.000000 

PROCESSOR  USAGE*  123.-?  UNITS 


Notice  that  H^(v;x,dn,T)  ■ 1 for  all  o .<  x _<  v + 8dn» 

The  Integral  (4.5)  Is  evaluated  numerically,  for  each  n ■ 2,3,...  on  a 
constant  grid  of  v values,  being  ^ ■ 1&,  1 ■ 0,  1,2,...,  M 
where  a ■ 25  [m/sec],  by  a formula  similar  to  (4.4).  In  Tables  2 and  3 
we  present  the  numerical  results  of  computing  five  exit  distributions 
recursively.  The  computations  were  performed  according  to  Program  BUL2 
given  in  Appendix  2.  The  difference  between  the  two  examples  Is  In  the 
value  of  the  a coefficient. 

5.  DETERMINATION  OF  THE  HIT/KILL  PROBABILITIES 

In  Section  2 we  Introduced  the  kill  probability  function  PK(vQ»R).  If 
the  bullet  goes  through  n trees  on  Its  way  to  the  target,  the  kill 
probability,  olven  the  last  exit  velocity  Vn  ■ vn  and  the  vector  of 
distances  , Is  PK(vn,  R-$n-T).  Accordingly,  the  kill  probability, 

given  N ■ n and  given  d^  Is,  for  trees  of  fixed  radius  T and  n >_  1 


*T) 


PK(x,R-?n-T)dHn(x;v0,d(n),T). 


(5.1) 


Furthermore,  since  the  conditional  distribution  of  the  points  of  location 

of  the  trees,  ■ Dj,  $2  * Oj+Dj**** »$n  " °l+***,Dn’ 
given  N ■ n , Is  like  that  of  ordered  statistics  from  a uniform  distri- 
bution on  (0,R),  the  conditional  kill  probability,  given  N ■ n and  T is 
for  n>l 
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TABLE  2.  Distributions  of  Exit  VI ocities;  vQ  3 1,300  [m/sec],  d2  3 200, 

d2  3 150,  d3  3 125,  d4  = 125,  d5  * 200  [m],  « 3 474,573.75  [m/sec2], 
g * .18,  T 3 .3  [ml. 


v/n 

1 

2 

3 

4 

5 

0. 

0.000000 

0. 000000 

O.TiOOOOO 

0. 000000 

0.000000 

29, 

0.000000 

0. 000000 

0. 000000 

0. 000000 

0.000000 

SO. 

0. 000000 

0. 000000 

0. 000000 

0. 000000 

0.000000 

79. 

0.000000 

0. 000000 

0. 000000 

0. 000000 

0.000000 

100. 

0.000000 

0.000000 

0. 000000 

0.000000 

0.000000 

123. 

0.000000 

0. 000000 

0. 000000 

0. 000000 

0.000000 

ISO. 

0.000000 

0. 000000 

0. 000000 

0. 000000 

0.000000 

ITS. 

0. 000000 

0. 000000 

0. 000000 

0. 000000 

0.000000 

200. 

0. 000000 

0. 000000 

0. 000000 

0. 000000 

0. 000000 

22S. 

0. 000000 

0. 000000 

0. 000000 

0. 000000 

0.000000 

2SO. 

0. 000000 

0. 000000 

0. 000000 

0.000000 

0.000000 

273. 

0. 000000 

0. 000000 

0. 000000 

0. 000000 

0.000000 

300. 

0. 000000 

0. 000000 

0. 000000 

0. 000000 

0.000000 

323. 

0.000000 

0. 000000 

0. 000000 

0. 000000 

0.000000 

330. 

0. 000000 

0. 000000 

0. 000000 

0. 000000 

0.000000 

373. 

0. 000000 

0. 000000 

0. 000000 

0.000000 

0.000000 

400. 

0. 000000 

0. 000000 

0. 000000 

0. 000000 

0.000000 

423. 

0. 000000 

0.000000 

0. 000000 

0. 000000 

0.003033 

430. 

0. 000000 

0. 000000 

0. 000000 

0. 000000 

0.009672 

473. 

0.000000 

0. 000000 

0. 000000 

0. 000000 

0.030402 

SOO. 

0.000000 

0.000000 

0. 000000 

0. 000000 

0.037226 

323. 

0. 000000 

0. 000000 

0. 000000 

0. 000000 

0. 082890 

330. 

0.000000 

0. 000000 

0. 000000 

0. 000000 

0.130717 

373. 

0.000000 

0.000000 

0.  000000 

0. 000000 

0.185662 

600. 

0. 000000 

0.000000 

0. 000000 

0. 000000 

0.244578 

623. 

0.000000 

0. 000000 

0. 000000 

0. 000000 

0.313746 

630. 

0.000000 

0. 000000 

0.  000000 

0. 000000 

0.413488 

673. 

0.000000 

0. 000000 

0. 000000 

0. 000000 

0.310973 

' 700. 

0. 000000 

0. 000000 

0.  000000 

0.021790 

0.599678 

723. 

0.000000 

0. 000000 

0. 000000 

0. 072218 

0.673214 

730. 

0. 000000 

0.  000000 

0.  000000 

0. 144083 

0.737767 

773. 

0. 000000 

0. 000000 

0.  000000 

0.223350 

0.823301 

800. 

0. 000000 

0. 000000 

0.  000000 

0.318064 

0.879620 

823. 

0. 000000 

0.  000000 

0.  000000 

0.430036 

0.920478 

830. 

0.000000 

0.  000000 

0.  000000 

0.380120 

0.949602 

873. 

0. 000000 

0.  000000 

0.  051706 

0.696111 

0.969486 

900. 

0.000000 

0. 000000 

0. 165134 

0.791303 

0.983231 

923. 

0. 000000 

0.  000000 

0.303905 

0.364773 

0.991339 

930. 

0. 000000 

0.000000 

0.443347 

0.916221 

0.993834 

973. 

0.000000 

0.  000000 

0.613173 

0.934470 

0.993143 

1000. 

0. 000000 

0.  000000 

0.734902 

0.977061 

0.999233 

1023. 

0.000000 

0.139957 

0.333303 

0.989343 

0.999713 

1030. 

0.000000 

0.394401 

0.925631 

0.995492 

0.999903 

1073. 

0. 000000 

0.603108 

0.96473 4 

0.998291 

0.999973 

1100. 

0.000000 

0.771304 

0.933291 

0.999433 

0.999993 

1123. 

0. 000000 

0.389233 

0.994716 

0.999840 

0.999999 

1130. 

0.236790 

0.933309 

0.998403 

0.  999963 

1.000000 

1173. 

0.647181 

0.984623 

0.999614 

0. 999994 

1.000000 

1200. 

0.832639 

0.993977 

0.999936 

1. 000000 

1.000000 

1223. 

0.940097 

0.999407 

1. 000000 

1.000000 

1.000000 

1230. 

0.992331 

1.000000 

l.  000000 

1. 000000 

1.000000 

1273. 

1. 000000 

1.000000 

1.000000 

1.000000 

l. 000000 

PROCESSOR  US POE:  20.3  UMITS 
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TASLE  3.  Distribution  of  Exit  Velocities;  vQ  = 1,300  [m/sec],  d1  * 200, 

d2  - 150,  d3  - 125,  d4  » 125,  dg  - 200  [m];  «-  1,000,000  Im/sec2], 
8 « .18,  T ■ .3  [mj. 


v/n 

1 

2 

3 

4 

5 

" ’o". 

o.  ooocnsTT 

0. 000000 

0.248346 

0.849444 

0.987763 

29. 

0. 000000 

0. 000000 

0.249499 

0.349883 

0.987313 

30. 

0.000000 

0. 000000 

0.232296 

0.831183 

0.987960 

79. 

0.000000 

0. 000000 

0.236781 

0.833297 

0.988273 

100. 

0.000000 

0. 000000 

0.262731 

0.836133 

0.988643 

129. 

0.000000 

0.000000 

0.269993 

0.839681 

0.989089 

130. 

0.000000 

0.000000 

0.273311 

0.863787 

0.989398 

173. 

0.000000 

0. 000000 

0.298032 

0.370212 

0.990161 

200. 

0.000000 

0. 000000 

0.313329 

0.376223 

0.990764 

223. 

0. 000000 

0.000000 

0.329014 

0.382340 

0.991492 

230. 

0. 000000 

0. 000000 

0.344104 

0.88837? 

0.992193 

273. 

0. 000000 

0. 000000 

0.371997 

0.396986 

0.992883 

300. 

0.000000 

0. 000000 

0.393368 

0.904323 

0.993623 

323. 

0. 000000 

0. 000000 

0.413887 

0.911426 

0.994328 

330. 

0.000000 

0. 000000 

0.448391 

0.920106 

0.994983 

373. 

0.000000 

0. 000000 

0.473916 

0.927413 

0.993663 

400. 

0. 000000 

0. 000000 

0.309298 

0.933481 

0.  996233 

423. 

0. 000000 

0. 000000 

0.333863 

0.942399 

0.996834 

430. 

0.000000 

0. 000000 

0.373010 

0.949899 

0.997338 

473. 

0.000000 

0.000000 

0.606748 

0.936373 

0.997804 

300. 

0. 000000 

0. 000000 

0.644192 

0.962792 

0.998202 

323. 

0.000000 

0.000000 

0. 676097 

0.968313 

0.998339 

330. 

0. 000000 

0. 000000 

0.714444 

0.973673 

0.998334 

373. 

0.000000 

0.000000 

0.744889 

0.973037 

0.999111 

600. 

0.000000 

0.000000 

0.781340 

0.982317 

0.999313 

623. 

0.  000000 

0.  022092 

0.813364 

0.983847 

0.999487 

630. 

0. 000000 

0. 037399 

0.340943 

0.988732 

0.999622 

673. 

0. 000000 

0. 132739 

0.868480 

0.991342 

0.999726 

700. 

0. 000000 

0.239789 

0.391901 

0.993404 

0.999806 

723. 

0. 000000 

0.310638 

0.912011 

0.993038 

0.999866 

730. 

0.000000 

0.392621 

0.930361 

0.996400 

0.999909 

773  • 

0.000000 

0.472333 

0.943940 

0.997427 

0.999940 

800. 

0. 000000 

0.343983 

0.933373 

0.998197 

0.999962 

823. 

0. 000000 

0.610063 

0.969236 

0.998779 

0.999976 

830. 

0. 000000 

0.679238 

0.977731 

0.999192 

0.999986 

873. 

0. 000000 

0.741204 

0.984292 

0.999481 

0.999992 

900. 

0.000000 

0.793373 

0.989221 

0.999677 

923. 

0. 000000 

0.843923 

0.992890 

0.999806 

0. 999998 

930. 

0. 000000 

0.337733 

0.993441 

0.999888 

0.999999 

973. 

0. 000000 

0. 922429 

0.997168 

0.999938 

om ****** 

1000. 

0.  087331 

0.947879 

0.998306 

0.999967 

1. oooooo 

1023. 

0.410648 

0.963621 

0.999033 

0.999984 

1.000000 

1030. 

0.364636 

0.978007 

0.999480 

0. 999993 

1.000000 

1073. 

0.676127 

0.986333 

0.999741 

0 # ww; 

1.000000 

1100. 

0.763202 

0.992334 

0.999883 

0. 

1.  oooooo 

1123. 

0.S32881 

0.993996 

0.999933 

1. oooooo 

1.000000 

1130. 

0.888612 

0.998139 

0 . 999984 

1. oooooo 

1.  oooooo 

1173. 

0.932262 

0.999302 

0 . 999996 

1. oooooo 

1 . oooooo 

1200. 

0. 964843 

0.999307 

0 _ 

1. oooooo 

1 . oooooo 

1223. 

0.986826 

0.999971 

i. oooooo 

1. oooooo 

1. oooooo 

1230. 

0. 998273 

1. 000000 

t.  oooooo 

1. oooooo 

l. oooooo 

1273. 

1. 000000 

1.000000 

1. oooooo 

1. OOOOOO 

1.000000 

PROCESSOR  US 
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where  the  vector  (^-T,  {g"  {fT $n-  ^ n-1  -2T ) Is  substituted  In  (5.1) 

for  d^.  For  n ■ 0 we  define  F.  (v_,  0,T)  * P|/(v,,,R).  Finally,  the  total 
kill  probability  Is  k 0 K 0 

0° 

FK(v0.T)  -e-R“i:-iSf-FK(Vn.T).  (5.3) 

We  discuss  here  two  numerical  procedures  for  the  estimation  of  fk(v0,T)  and 
their  accuracy.  The  values  of  the  function  PK(x,  R-?n-T)  are  given  or 
computed  on  the  same  grid  of  points  * 1a,  1 * o,...,25,  as  that  used 

for  the  numerical  computation  of  the  functions  Hn(v;vQ,d^n\T).  Thus, 
the  function  (5.1)  Is  evaluated  numerically  according  to  the  formula 

♦K(Vo,n,d^,T)»  (5.4) 

m 

£ PK(«7?,R-VT)CHn^1*vo4(n)*T)  " Hntol-l;V4(n)*T)]'  n^* 
1-1 


It  Is  much  more  complicated  to  evaluate  the  function  FK(v0.n,T)  numerically. 
We  have  Introduced  here  two  methods  for  a Monte  Carlo  estimation  of  (5.2). 


Method  I 


For  each  n,  n ■ 1,  2,...,  NP,  perform  the  following  Monte  Carlo  estima- 
tion Independently.  Simulate  n uniform  random  variables  on  (0,R).  Let 

be  the  order  statistics  of  these  simulated  variable.  Compute 

the  function 'frK(v(),n,d^n^ ,T)  for  these  values.  Repeat  this  simulation 
Independently  NS  times  and  average  the  resulting  ♦k(*)  values.  This 

average  £K(vQ,n,T)  Is  an  unbiased  estimator  of  ^(vQ,n,T).  Let  be  the 
sample  variance  of  the  simulated  ^k(*)  values.  An  estimator  of  the  variance 
of  fyv0,n,T)  Is  thus  S^/NS.  Finally,  the  kill  probability 

FK(vo.T)  is  estimated  by 


Pk<Vt>  ’ 


NP  / D\n  a 

£ ^n^~  FK(Vn,T** 

n«0  n*  * 0 


(5.5) 


fyvo.Q.T)  - Fk(vo.0,T)  - Pk(v0.R). 


NP  Is  a sufficiently  large  Integer,  so  that  the  sum  of  the  Poisson  pro- 
babilities, for  n larger  than  NP,  Is  sufficiently  small.  In  the  following 
numerical  examples  we  specified  the  value  of  NP  - INT  (uR+3\/uR),  where 
INT(x)  denotes  the  Integer  smaller  than  or  equal  to  x.  Let  pos(n;  uR) 
denote  the  Poisson  probability  of  N ■ n,  with  parameter  uR.  The  variance 
of  ^«(v0»T)  Is  estimated  by 


The  square  root  of  (5.6)  Is  the  standard-error  of  the  estimate  of  the  kill 
probability. 


Method  II 

According  to  the  second  method  Me  estimate  fK(vQ,T)  In  one  simulation 
string;  repeat  this  estimation  Independently  NS  times  and  take  the 
average  of  the  Individual  estimates.  More  specifically.  Starting  with 

fk(v0,°,t)  ■ pic<v0»R)  set 


*lS1)(VT>~ * fk(v0,o,t). 


We  then  set  n ■ 1 and  generate  a uniform  random  variable  from  (0*R).  This 
value  Is  set  equal  to  d1  and  the  function  ^(v^l.d^.T)  Is  computed. 

The  value  of  ^^(V0»T)  Is  then  changed  to  $^^(v0,T)  + pos(l;Ru)* 
♦K(v0,l,di,T).  We  set  then  n ■ 2,  generate  a uniform  random  variable  from 
(O.R-dj),  which  Is  set  to  be  equal  to  dg.  We  then  compute  ♦K(vQ, 2, d^dg.T) 

and  set  ^k^^(v0,T)— Pk^^(v0,T)  + pos(2;MR)* '*'K(v0,2,d1,d2,T).  This 
simulation  algorithm  Is  continued  until  n ■ NP.  According  to  this 
algorithm 


i designates  the  index  of  the  simulation  run,  i = 1,...,NS.  Finally, 


Fk(Vq,T)  Is  estimated  by 


» l*  f,  ^(1)(vft,T). 


NS 


(5.8) 


An  estimate  of  the  variance  of  ”„(vrt,T)  Is  obtained  by  computing  the 

sample  variance  of  the  ^^(V0»T)  values  (1  ■ 1,...,NS)  and  dividing 

this  sample  variance  by  NS.  The  square-root  of  this  variance  Is  the 

standard  error  of  ^(v0,T).  In  Appendix  3 we  provide  Program  BUL6  which 

estimates  the  kill  probabilities  ^K(v0,T)  according  to  Method  I.  Program 

BUL7  given  in  Appendix  4 is  designed  for  the  estimation  of  F^0J) 
according  to  Method  II.  In  both  programs  we  considered  the  special 
function 


' t 


if  v0-$R  > vk 


otherwise. 


(5.9) 


In  Program  BUL6  this  special  P^( vQ ,R ) function  Is  programmed  In  the  main 
routine.  Program  BUL7  Is  written  In  a more  general  manner.  The  function 
Pk(Vq,R)  is  computed  as  a subroutine  function  and  can  be  changed  without 
altering  the  main  program.  In  Table  4 we  provide  estimates  of  the  kill 
probability  computed  according  to  Method  I and  Method  II  with  the  function 
Pk(v0»R)  given  by  (5.9).  For  Method  I we  present  the  average  estimates 

nf  \b  ( u n r|(n)  Tl  an/t  tka  rnrracnnnHinn  omnia  etanHarH  s4a\H  Atlnnc  rrf 


of  *|<(v0»n,4'  »T)  and  the  corresponding  sample  standard  deviations  of 

these  estimates  for  n ■ 0,1,...,NP.  The  estimate  ^( V0»T)  and  Its 
standard  error  are  given,  too.  For  Method  II  we  present  the  estimates 


(VT>  ^or  * * average  fy(vQ,T)  and  Its  standard 

error.  We  see  that  the  accuracy  of  Method  I Is  somewhat  higher  than 
that  of  Method  II.  However,  Method  I requires  more  than  5 times  longer 
computer  time  than  Method  II.  It  seems  justified  to  recommend  the  use 
of  Method  II. 


Ik 


1 
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TABLE  4.  Estimates  of  the  Kill  Probability  by 

Method  I and  Method  II;  vQ  ■ 1,300  Im/secJ, 

u • .005,  R ■ 1,000  [m],  a • 474,573.75  [m/sec2], 

8 ■ .18,  ■ 500  Im/secJ. 


METHOD  I 

METHOD  II 

n 

♦l^(v0,n,lt^n^  *T) 

s„ 

1 

#K°  <V  T) 

0 

1.000000 

0.000000 

1 

0.848318 

1 

1.000000 

0. 000000 

2 

0.832843 

2 

1.000000 

0.000000 

3 

0.833372 

3 

1. 000000 

0. 000000 

4 

0.831003 

4 

0.888848 

0. 001109 

3 

0.813378 

3 

0.788708 

0. 092821 

8 

0.813843 

8 

0.293843 

0. 039027 

7 

0.804393 

? 

0. 082802 

0.019183 

3 

0.831333 

8 

0.015848 

0. 003773 

9 

0.884387 

9 

0.001833 

0. 000373 

10 

0.892343 

10 

0.000413 

0. 000234 

11 

0.000033 

0. 000009 

A 

'I* 

.633280 

PK  » .638706 

S.E.lfy  ■ .005485 

S.E.{fK> 

» .008206 

Processor  Usage:  1977.9  units 

Processor  Usage:  368.5  units 
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6.  GENERALIZATION  FOR  RANDOM  TREE  RADIUS 


In  the  present  section  we  Indicate  how  the  previous  results  can  be 
generalized  to  the  case  of  randan  trunk  radius,  T.  Let  F^(T)  be  the 

c.d.f.  of  T.  Given  that  N » n,  the  corresponding  values  Tj Tn  are 

Independent  and  Identically  distributed.  The  following  modifications 

are  needed.  Let  T^  » (T^ Tn).  The  exit  velocity  distributions 

are  then  computed  as 


Hl(v*vo*dl,Ti)  , for  n * 1 (6.1) 

H2(v;Vdl*d2,Tl*V  *y*Hi(v;x,d2*VdHl(x;VdlTl)  . for  n * 2 (6.2) 


and 


H (vv  d(n)  T(n);  */Hi(vix*dn»Tn)dHn-l(x;vo»^(n”1)»I(n“1))*  for  n > 2*  (6*3) 

nv  ’ o*~  J 


Given  these  exit  velocity  distributions,  we  determine  for  each  n > 1,  the 
conditional  kill  probabilities 


*K(v0,n,d("\l(n))  * /PK(x.R-«)dHn(x’vo *i(n)*I(n))* 

The  kill  probabilities,  given  vQ,  and  |N*n|  are  computed  then  as 

Vvo*n)  -/-/WW  dW» 


(6.4) 


(6.5) 


where  Pk^o'0^  is  a 9enera1  ization  of  the  function  defined  in  (5.2) 

obtained  by  Integrating  *K(v0,n,  ^-T,  f2-  ^-21 *n-  jn_r2T, 

over  the  simplex  of  Jj_<  {.><.. • • ._<{n.  This  Integral  Is  estimated,  as  In  the 
previous  section  by  Method  I.  The  Integral  (6.5)  can  then  be  evaluated 
numerically  either  by  evaluating  a discrete  version  of  It,  or  by  simula- 
tion. If  one  employs  Method  II  then  (6.5)  is  evaluated  for  each  n In  the 
same  string  of  computations.  Finally,  in  Method  I the  functions  F^Q) 

are  computed  as  Poisson  averages  of  F,Jv  .n). 

^ 0 
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APPENDIX  I.  FORTRAN  PROGRAM  BULLET 
(Computation  of  H2(v;  vQ,  d^,  d2>  T) 


100 

103 

no 

120 

130 

130 

1*0 

170 

iao 

190 

200 

210 

220 

230 

243 

24« 

247 

230 

260 

270 

280 

281 

282 

283 

284 
283 
286 
237 
288 

289 

290 

291 
293 
300 
310 
320 
330 
340 
330 
360 
370 
330 
390 
400 
410 
420 
423 
430 
440 
430 
460 
470 


DIMENSION  M<3> 
V0-1300. 
01*200. 

02*230. 

A— 1730000. 
>-.18 
Vl-V0->*01 
V2*Vl->*02 
N-300 


VO  4-  v 

0 

A 4.  9 
B 4.  g 


K*13 

DO  1 1*1  »K 
AI-I-1 
VI-100. *AI 
00  20  L-1.3 

T 4-  trunk  radius 

i ® • 1 •HL 

Wl-l . - < (V2-V2-VI *VI>  ' <2.  *A*T)  > **2 
IF<WI>  2,3.3 

2 Wl-O. 

3 H<D-SQRT<UI> 

WlS«Vl*Vl-2.*A*r 

if<wis>  40>4i.4i 

40  W1S-0. 

41  VIS-SQRT(WIS) 

VS*V1S->*D2 

IF <VI-VS)  42.43,43 

42  RI-V2-VS 
YS-VS 

SO  TO  30 

43  R I *V2-V I 
YS-VI 

30  00  4 J-l,H 
AJ-J 

YJJ-Y  J— R I/'AM 
YJH* <YJ*YJJ> ^2. 

ZJ-  < <YJ*YJ-VI*Vp C2.  ♦4-T1  > **2 
ZJJ-  < <Y  J J4Y J J-V I • V I > /■  <2.  *A*T)  1 **2 
IF  <ZJ-1.)  3.4,4 

3 IFCJJ-U.)  6,4.4 

6 S J-SQRT  < 1 . -Z  J> -SORT  < 1 . -2 J J) 

UJ»<  CVI*V1-<YJM*B*02>  **2>  s (.2.  *A*T>  ) **2 
IF  OJJ-1.)  7,4,4 

7 T J-SQRT <1.-UJ> 

H(U)-H<U-TJ*6J 

4 CONTINUE 
20  CONTINUE 

PRINT  10.VI, <H<L>.C*1,3> 

10  FORMAT <3X. F6. 0»3F 10.6) 
l CONTINUE 
• END 


< 


i 


APPENDIX  II.  FORTRAN  PROGRAM  BUL2 
(Recursive  computation  of  Hn(v;  vQ,  d^,  T) 


OOOIO 

DIMENSION  H<3*52>  *D<3>  *U<3>  , 

00020 

DATA  <D  <I>  » I»1 » 5) '200. * 150. * 125. . 125. » 200. ✓ 

00030 

2-25. 

00040 

K*3 

00030 

V 0*1300. 

00060 

3*1000000. 

00070 

T*.3 

00080 

B*.  18 

00090 

M-52 

00100 

flH*H 

00110 

vi-vo-o<i)-i 

00120 

U<1>*V1 

00130 

do  ioo  i*i*m  , 

00140 

AI-I-l 

00150 

VX*AI«£ 

00160 

21*  <V1*V1-VI*VI5 / <2. *A*T5 

00170 

IF  <2I>  101.101*102 

00130 

101 

H<1  * l>  *t . 

00190 

60  TO  100 

00200 

102 

IF <21-1. > 103.103*104 

00210 

103 

H<1»  D-SQRT <1. -21-21) 

00220 

60  TO  100 

00230 

104 

H<1»  I)  *0. 

00240 

100 

CONTINUE 

00250 

DO  200  L-2.K 

00260 

U <U  «U  0.-0  -B  <U  *8 

00270 

VL-U<U 

00230 

DO  300  I*1*M 

00290 

AI-I-1 

00300 

VI-AI-S 

00310 

IF<VI-VU  301*301*302 

00320 

301 

H<L»  I)  *0. 

00330 

do  400  j*i*n 

00340 

AJ-J-1 

00350 

XJ-£*<AJ*.3> 

00360 

2J*  <XJ-D <U -B> ♦ <XJ-D  <U  *B> -VI-VI 

00370 

ZJ"ZJs  <2. 

00330 

IF  <2J>  401*402*402 

00390 

401 

PIJ-1. 

00400 

60  TO  403 

00410 

402 

IF  <2J-1.>  403*403*404 

00420 

403 

PIJ-SQRT  <1.-2J-2J> 

00430 

60  TO  403 

00440 

404 

PIJ*0. 

00430 

405 

IF <J-1>  406*406.407 

00460 

406 

H <L.  I>  -H  <L*  I>  *P I J-H  <L-1 . J> 

00470 

60  TO  400 

00430 

407 

H<L.  I)  -H  <t*  I>  *PI J*  <H  <L-1  * J>  -H  <L-1  * J-l>  ) 

00490 

400 

CONTINUE 

00500 

60  TO  300 

00510 

302 

H<L»  I)*l. 

00520 

300 

CONTINUE 

00S30 

00540 

00550 

00560 

00570 

00530 

200 

CONTINUE 

DO  500  I — 1 * M 

ix-x-l  .c 

ai-ii  , c ^ 

VI-AI-E  TH^S  ^ fj 

PRINT  510.VI. <H<L* I> *L«l.tO  f&P  ^ 

00590 

510 

FORMAT <3X>F6.0*3F10.6> 

00600 

500 

CONTINUE 

00610 

END 

E «■  A 
A a 
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APPENDIX  III.  FORTRAN  PROGRAM  BUL6 
(Estimation  of  f ^(.vQ , T)  by  Method  I) 


ioosuii,ramdx> . • 

no  DIMtMSIOM  <3v31> -Hv31>  »0<100>  >U<100> 

110  ft-1000. 

130  V 0*1300. 

140  3*474973. ?3 

190  AU-3. 

140  AMU-R/AL 

170  T-.3 

110  VK-300. 

190  B*. 13 

100  MU-100 

110  MF*ll 

no  n-si 

130  AH-M 

140  MS- 10 

ISO  AMS -MS 

140  K-13. 

170  Y-AAMB<-1.> 

110  00  1 t-l.MR 

190  Y-RAMB<0.> 

300  1 COMTIMUt 

310  K-0 

310  FK-AOSOt.AU 

330  VR-VO-I-A 

340  IF  <VR-VtO  11 >11 >12 

390  11  FK-O. 

340  «l  TO  13 

370  11  FK-1. 

330  13  TFK-FK-FK 

390  3FK-0. 

400  FRIMT  14>K>FK..QFK 

410  14  FORMAT <SX> 14. 2F10.4) 

410  00  400  K-l.MF 

430  SFK-O. 

440  SSFK-O. 

490  DO  30  US-1 > MS 

440  U-0. 

470  Q-0. 

430  00  3 J-l.K 

490  Y-RAMD  v 0.  > 

900  0<J>-Y-«.R-W» 

310  W-V-D<J> 

310  3-3*B*0<J) 

330  UU>-4 


•[library  function  for 
generating  uniform  M.V.'s 
J>n  (0,1) 


A 4-  a 
AL  4.  RU 


This  loop  Is  just  for 
generating  the  first 
NR  • 100  uniform  R.V. 's 
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APPENDIX  III  (Continued) 

340 

i coMTimji 

330 

IF<K-l>  3.3.9 

340 

3 VI«V\MJ<.1> 

370 

VKK«VK-*«M-U<K> 

330 

JX«<VI*V1-VKK*VKK>/ <3. *4*0 

340 

lF<BO  31.33.33 

400 

31  FK-O. 

410 

30  TO  33 

MO 

33  XF<»-1.>  33.33.34 

430 

33  PK«l.-$QFT <!.-»♦»> 

440 

30  TO  93 

430 

34  PK«t. 

440 

33  SFK-3FKMK 

470 

SSFK-SSFKMKMK 

MO 

'50  TO  30 

490 

9 V1«V0-<JU> 

700 

DO  too  t«i»n 

710 

*X«t-l 

730 

v i-« -Mr 

no 

2i»tvi*vi-vi«vi>>'<i.«**n 

740 

if<si>  101. 101.103 

730 

tot  9 <!>■!. 

740 

30  TO  100 

770 

103  IF<ZX-1.> 103* 103. 104 

730 

103  3 <X> -SORT <1.-31 *31) 

790 

30  TO  100 

300 

104  3<:>>o. 

• 10 

100  CONTINUI 

330 

SO  300  L«3»K 

330 

VU«V0HJ<U 

340 

DO  300  I«l.n 

330 

*X«I-l 

340 

VX-CMI 

370 

IF<VI-VO  301. 301. 303 

330 

301  H<1>  »0. 

390 

00  400  J-l.M 

900 

910 

930 

2J" <XJ-D  <U> •» • tX J-D vL>  *1) -V l *V I 

930 

3J"3J'<3.M*r» 

940 

IF <ZJ>  401.403.403 

930 

401  F1J-1. 

940 

30  TO  403 

970 

403  IF <2J-1 .>  403.403*404 

1 

990 

403  FtJ«S4FT<l.-2J*3J> 

990 

30  TO  403 

1000 

404  FU-0. 

toto 

403  IF < J-l>  404.404.407 

1030 

404  M<X>M<l>M!J*3vJ> 

1030 

30  TO  400 

1040 

407  H<l>«H<tNMlJ*<3vJ>-3vJ-n> 

1030 

400  COMTINU* 

1040 

10  TO  300 

1070 

503 

1030 

300  COMT1HUI 

1090 

00  300  l«l.H 

1100 

3 <!>•«*!> 

1110 

1130 

300  COMTlfH.'i  

300  CCMTIMJf  THIS  fX«  Ui  QUAI.I “ r*01™ 

mm 

J»VY  ntHJH  ' S'ifl  ft' DDC 

1 

1 

III-2 
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HP***—-"  •»  f£3 

APPENOIX  III  (Continued) 


1130 

1140 

1130 

11*0 

1170 

1100 

11*0 

1*00 

1*10 

12*0 

1*30 

1*40 

1*30 

12*0 

1*70 

1*10 

12*0 

^300 

1310 

13*0 

1330 

1340 

1330 

13*0 

1370 

13*0 

13*0 

1400 

1410 

14*0 

1430 

1440 

1430 

14*0 

1470 

14*0 

14*0 

1300 

1310 


1340 

1330 

13*0 

1370 

13*0 

13*0 

1*00 

1*10 

1*20 

1*30 

1*40 

1*30 

1*40 

1*70 

1M« 

1**0 

1700 

1710 

17*0 


VKK-VK**-4-U<K> 
IK-INT<VKK«*>fl 
4K-1.-<H<XK>-H<XK-1>>/*. 
S4K-S4K44K 
SS4K-SS4K*4K«4K 
30  COhTINUt 
44K-34K^4H3 

S8K-<ftnS*$S4K-$4K-$4K> '4HS 
S4K«Sa*T  (SBKs  <4tlS-l . i > 

44 1 NT  34.K.44K. S4K 
4K-4QS  <K»*0 -4QS  (K-t.AO 
T4K-T4K  44K-44K 
84K-i84K'*4K*4K-S4K*$4K^4NS 
400  CONTimil 

SB»SQ4T<Q4K> 

44XHT  4l.T4K.SB 


41  4XBW4T<-".3X.rH«S.  4K«.410.4»7HCS. 

4UNCTtOM  403<J»*O  _ 

r-j 

I4<».ij*.  10.>  SO  TO  9 
14  <1>  1.3.3 
1 403-0. 

00  TQ  10 
3 4Q3-CX4  <-l) 

SO  TO  10 

3 40S-«X4<-1> 

4-403 

BO  4 K-l.I 
4K-K 

4-4-1X4K 
403 -403 >4 

4 CCHTXMU* 

SO  TO  10 

9 4I-I-.3 

ZI-<«I-»>'SQ*T<*> 

40S-CNBX  <ZX> 

10  RITURM 

m 

4UMCTX0M  CMBX  OO 
Y-X 

ISUTCH-0 
14  <Y>  1.3.3 
1 Y-41S  <Y) 

ISUTCH-1 
3 4-. 331441* 

. I1-. 31*39133 
**—  .33434379 
*3-1.7914-7* 

14—1.931233* 

•3-1.3302744 

T-l.y<l.-4-Y) 

»-. 3*8*423 *€.X4 OY-Y'2. ) 


S8-*410.«> 

Subroutine  function  for 
computing  the  Poisson 
c.d.f. 


For  large  AL  it  needs 
the  normal  c.d.f. 


Subroutine  function  for 
computing  the  standard 
normal  c.d.f. 


OMBX-l. -4*<*l*T4*3-T4r*I3-r-T-T*t4-^T**4»4*5-<T*-3) > 


3 

4 


X4<XSUTCM>  3.4.3 

QMBX-l.-*MBX 

cnox-ariDx 

ornjwi 

mo 
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APPENDIX  IV.  FORTRAN  PROGRAM  BUL7 
(Computing  ^CvQ»  T)  according  to  Method  II) 


UttPLSa’,"ffltaiian  e<3*>  .h<3«  ,o<i 

00030  Ml  000. 

00040 

11*. 0*3 

00090 

V0-1300. 

000*0 

*•474973. 79 

00070 

AL-H.-0 

000*0 

T-.13 

00090 

VK-300. 

00100 

1* 

00110 

NMlOO 

ooiao 

HMIMT  (W.O.  *SQ7T  «HJ  > 

00130 

n-3* 

00140 

00130 

NS»10 

001*0 

W1S-MS 

00170 

*•23. 

001*0 

y»**nd<-i.> 

00190 

oa  i i*i,nr 

00*00 

Y-**ND<0.> 

00*10 

1 CONTI MU* 

00*20 

K»0 

00*30 

7K«7K<V0«*> 

00*40 

T7K«*K**0S <K,*L> 

00*90 

U7K-T7K 

00*1,0 

Q7M0. 

00*70 

S7K-0. 

00*00 

00  20  US* 1 , NS 

00*90 

T7K-W7K  • 

00300 

U-0. 

00310 

0*0. 

003*0 

Ml 

00330 

Y«#*NO  <0. ' 

00340 

0<1>»Y* <*-40 

00330 

W"W*D<1> 

003*0 

Q*G*t*D<l> 

00370 

U<l>-« 

003*0 

V1»V0-U<1> 

00390 

4K-0. 

00400 

do  ioo  i»i»n 

00410 

*I«l-i 

004*0 

vt-c«*t 

00430 

2><vt*vi-vi*vi>/<a.***T> 

00440 

17  <21>  101,101,10* 

00430 

101  *<!>•!. 

004*0 

80  TO  103 

00470 

10*  17  <21-1 . > 1 03, 103* 1 04 

004*0 

103  8<1>«SQ*T <1. -21*21) 

00490 

80  TQ  109 

00900 

104  8<!>aO. 

00910 

103  VlT«Vl*€''». 

009*0 

17<1.3T.1>  80  TO  106 

00930 

7K»*K*7K  <VIT,*-U1  *6il> 

00940 

80  TO  100 

Library  function  for 
generating  uniform 
R.V.'s  on  (0,1). 


BL  - u 
AL  «■  uR 
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APPENDIX  IV  (Continued) 


00330 

104 

<VIT»  F-U>  * <S  <I>  > 

00340 

too 

CONTINUE 

00370 

TFK-TFK-FK-.FOS  <K.FL.' -FOS  <K-1 . FO > 

003*0 

00  300  L-3.NF 

00390 

Y-FFND<0.> 

00400 

D <u  *y+  <F-U> 

00410 

W-W-D<L> 

004*0 

Q"G+B*D <L> 

00430 

U<L>-0 

00440 

VL-VOHJ  <U 

00430 

PK-0. 

00440 

DO  300  1*1 »M 

00470 

FI-l-1 

004*0 

VI -1*4 I 

00490 

VIT-VI*«'3. 

00700 

IF<VI-VO  301.301.303 

00710 

301 

H<I>-0. 

007*0 

DO  400  J-l.N 

00730 

FJ-J-1 

00740 

XJ««*<4J*.3> 

00730 

2 J-  <X  J-D  <L> ♦*> ♦ vXJ-D <U ♦*) -V I *V I 

00740 

2J-ZJ'  <3.  **+T> 

00770 

IF  <2J>  401.408.403 

00790 

401 

FIJ-l. 

00790 

SO  TO  403 

00900 

40* 

IF<2J-1.)  403.403.404 

00*10 

403 

F X J-SQRT  < 1 . -2 J*2 J) 

oo*to 

SO  TO  403 

00*30 

404 

FIJ-O. 

00*40 

403 

IF CJ-1)  404.404.407 

00*30 

404 

H<I>-H<I>  -FI J-S<J> 

00*40 

SO  TO  400 

00*70 

407 

M<I>»M<I>*FIJ*<S<J>-S<J-1>> 

00890 

400 

CONTINUE 

00990 

80  TO  303 

00900 

303 

H<X>-1. 

00910 

303 

FK-FK <VIT.F-W> 

00980 

304 

IF<I.ST.1>  SO  TO  304 

00930 

00940 

SO  TO  300 

00930 

304 

PK«^ICH»K*iH<I)-Ha-l>) 

00940 

300 

CONTINUE 

00970 

DO  300  1*1. M ' 

00990 

9<X>-H<I> 

00990 

300 

CONTINUE 

01000 

tfk-tfk-fk* <fos  <L»  FO -FOS  <L-1 . FO > 

01010 

*00 

CONTINUE 

01080 

sfk-sfk*tfk 

01030 

qfk«qf*otpk*tfk 

01040 

FFINT  S4.LS.TFK 

01030 

34 

FOFNFT <3X. I4.F10.4) 

01040 

30 

CONTINUE 

01070 

FFK«SFK^FN$ 

01090 

VFK"  <FNS-GFK-SFK*SFK>  ✓FNS 

01090 

SDK-SOFT <VFK'  vFNS-l. ) 

01100 

FFINT  41.FFK.SDK 

OHIO 

41 

FOFNFT  <3X.F10.6.3X<F10.4) 

011*0 

END 

APPENDIX  IV  (Continued) 


01130 

7VIMCTIQM  wreu.ao 

Subroutine  function  for 

01140 

01130 

I«J 

computing  the  Poisson  c 

011*0 

17  <t.  52.  10.)  SO  TO  3 

01170 

I7«>  1.2.3 

011*0 

i tos-o. 

011*0 

80  TO  10 

ouoo 

2 7<K-«X7<-0> 

01*10 

80  TO  10 

01220 

3 70S»«X7<-» 

01230 

7-7QS 

01240 

DO  4 K-t»l 

01230 

mm* 

012*0 

01270 

70S-7GS*7 

012*0 

4 coriTimic 

012*0 

80  TO  10 

01300 

3 *!•!♦. 3 

01310 

ZI« <«!-»> ^SQ*T(I> 

01320 

70s«cndx<2j> 

01330 

10  WTU7M 

01340 

040 

01330 

FUNCTION  CN0X(X> 

Subroutine  function  for 

013*0 

01370 

Y«X 

tSUTO4*0 

computing  the  standard 

013*0 

IP  CO  1.2.2 

normal  c.d.f. 

013*0 

1 Y«**S(Y> 

01400 

tSUTCH*! 

01410 

2 pm. 231*41* 

01420 

01*.  31)33133 

01430 

12***  • 33*3*373 

01440 

*3«1. 7*1477* 

01430 

14—1.3212339 

014*0 

03-1.33 02744 

01470 

t«i./<i.4*4y> 

0t4*0 

Pm.  3989423«<X7  <-Y*Y^2. ) 

014*0 

OMOX-l .**4<*1*T*0*4T4T*I3*T4T4T*044<T4*4> ♦l3*<T—3>) 

01300 

I7<ISUTCH>  3.4.3 

01510 

3 ONOX-t.-ONOX 

01320 

4 CMOXwJHOX 

01330 

Rrruon 

01340 

040 

Subroutine  function  for 

01330 

7UMCTXW4  7K(V.7> 

013*0 

01370 

u«v 

Pm* 

computing  PK(vQ»  R) 

013*0 

(*(■300. 

013*0 

*•.13 

01*00 

U$«M-7*» 

01410 

17  (MS  “4*0  1.1.2 

01*20 

t 7K-0. 

01*30 

SO  TO  10 

0t*40 

3 TK-t. 

01430 

10  72TU7N 

014*0 

040 

Will  PAOt  LS  B&SI  QULUUTY  HMgXK 
am  «rt  MWSUK  SO  BOO 
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